import os
import git
# Our numerical workhorses
import numpy as np
import scipy as sp
import pandas as pd
# Import Interactive plot libraries
import bokeh.plotting
import bokeh.layouts
from bokeh.themes import Theme
import holoviews as hv
import hvplot
import hvplot.pandas
import panel as pn
# Import the project utils
import covid19
bokeh.io.output_notebook()
hv.extension('bokeh')
# Set PBoC style for plot
theme = Theme(json=covid19.viz.pboc_style_bokeh())
bokeh.io.curdoc().theme = theme
hv.renderer('bokeh').theme = theme
In this notebook we'll explore a simple SIR-like model to better understand the effect that an increasing testing capacity would have in isolating the COVID-19 cases. Here we'll present interactive plots in which parameters can be modified via sliders. This interactive feature facilitates the exploration of the complicated parameter space.
For this model we will use the well established SIR-like model in which the population is compartmentalized into the following groups:
Note that we do not distinguish between people that recovered form the disease or that passed away. As far as the model goes, there is no distinction to how the dynamics of the epidemic change. Nevertheless if needed the model could easily divide the $R$ population into these two groups.
The assumption of our model are the same as the typical SIR model. What this means is that we assume we operate at the limit of large populations where small fluctuations are averaged out, so we use deterministic ODEs. We also assume a well-mixed population where contacts follow mass-action dynamics. All these assumptions can be easily relaxed, but for now we just want to get the general possible scenarios that can result from running the dynamics with different parameter regimes.
Let's now write down the equations for each of the populations. First for the susceptible people we have
$$ {dS \over dt} = - \beta_S S I_S^{(q)} - \beta_A S I_A - \beta_A^{(q)} I_A^{(q)}, $$where $\beta_s$, $\beta_A$, and $\beta_A^{(q)}$ are the infective rates for each of the infective populations.
$$ {dE \over dt} = \beta_S S I_S^{(q)} + \beta_A S I_A + \beta_A^{(q)} I_A^{(q)} - \gamma_S E - \gamma_A E, $$where $\gamma_S$ and $\gamma_A$ are the rates at which exposed people become infective presenting symptoms or being asymptomatic, respectively. For the three infected populations we have
$$ {dI_S^{(q)} \over dt} = \gamma_S E - \mu I_S, $$$$ {dI_A \over dt} = \gamma_S E - \theta I_A - \mu I_S, $$$$ {dI_A^{(q)} \over dt} = \theta I_A - \mu I_S, $$where $\theta$ is the rate in which an asymptomatic person goes into quarantine. This rate will depend on the availability of tests, the accuracy of such tests, and the probability that a person actually complies with the quarantine. We also assume that all infected groups either recover or die at the same rate $\mu$. With this the dynamics for the removed group are of the form
$$ {dR \over dt} = \mu \left(I_S^{(q)} + I_A + I_A^{(q)} \right) $$To look at the dynamics we will numerically integrate the system of ODEs using scipy's ode_int function. For this we first need to define a function that computes the right-hand side (rhs) of
def system_rhs(x, t, beta_s, beta_a, beta_aq, gamma_s, gamma_a, theta, mu):
"""
Right hand side of system. Returns time derivatives
Parameters
----------
x : array-like. len(x) = 6
Array containing all 6 compartmentalized populations in the following
order: S, E, I_s, I_a, I_aq, R.
t : array-like.
Time array in which to evaluate the dynamics.
beta_s, beta_a, beta_aq : float.
Infection rates for each of the three infected populations.
gamma_s, gamma_a : float.
Rate in which exposed individuals become infective with symptoms or
asymptomatic, respectively.
theta : float.
Rate in which an asymptomatic individual goes into quarantine.
mu : float.
Rate in which infected individuals are removed from this group.
Returns
-------
system_rhs : array-like.
Returns evaluation of the right-hand side of the ODE system
"""
# Unpack different population groups
S, E, I_s, I_a, I_aq, R = x
# Compute time derivatives
dS_dt = -beta_s * S * I_s - beta_a * S * I_a - beta_aq * S * I_aq
dE_dt = (
beta_s * S * I_s
+ beta_a * S * I_a
+ beta_aq * S * I_aq
- gamma_s * E
- gamma_a * E
)
dI_s_dt = gamma_s * E - mu * I_s
dI_a_dt = gamma_a * E - theta * I_a - mu * I_a
dI_aq_dt = theta * I_a - mu * I_aq
dR_dt = mu * (I_s + I_a + I_aq)
# Return results as numpy array
return np.array([dS_dt, dE_dt, dI_s_dt, dI_a_dt, dI_aq_dt, dR_dt])
To test that we can get the system up and running let's choose some simple initial condition and random set of parameter values. We'll then numerically integrate the system using scipy.integrate.odeint().
# Define time (in a.u.) for integration
t = np.linspace(0, 15, 500)
# Set initial conditions for population
init_infection = 1e-2
x0 = np.array(
[
1 - init_infection, # S
0, # E
init_infection / 3, # I_s
init_infection / 3, # I_a
init_infection / 3, # I_aq
0, # R
]
)
# Arbitrary set of parameters
beta_s = 1
beta_a = 5
beta_aq = 1
gamma_s = 1
gamma_a = 3
theta = 1
mu = 0.5
# Pack parameters into tuple
args = (
beta_s,
beta_a,
beta_aq,
gamma_s,
gamma_a,
theta,
mu,
)
# Numerically integrate system
x = sp.integrate.odeint(system_rhs, x0, t, args=args)
Let's now show the solution with an interactive Bokeh plot.
# Set colors for curves
colors = bokeh.palettes.d3["Category10"][5]
# Split populations
S, E, I_s, I_a, I_aq, R = x.transpose()
# Setup plot
p = bokeh.plotting.figure(
plot_width=500,
plot_height=300,
x_axis_label="time (a.u.)",
y_axis_label="fraction of population",
)
# Populate plot with glyphs
p.line(t, S, line_width=2, color=colors[0], legend_label="S")
p.line(t, E, line_width=2, color=colors[1], legend_label="E")
p.line(t, I_s, line_width=2, color=colors[2], legend_label="I_s")
p.line(t, I_a, line_width=2, color=colors[3], legend_label="I_a")
p.line(t, I_aq, line_width=2, color=colors[4], legend_label="I_aq")
p.line(t, R, line_width=2, color=colors[4], legend_label="R")
# Place legend
p.legend.location = "center_right"
# Show plot
bokeh.io.show(p)
The results follow the general SIR-like dynamics in which infected people initially increase exponentially and then decrease as the number of removed individual increases.
With a large parameter space let's turn our attention into making a set of interactive plots to explore this space.
panel¶Now that we know that the system is behaving as we would intuitively expect, let's set the parameters such that we can control them via panel. This is a new library that allows easy generation of interactive widgets.
To have a better handling of the model in terms of intuitive parameters we will split some of the parameters into more sub-parameters that we can freely explore. For example the rate $\theta$ at which asymptomatic infected individuals go into quarantine will be written as a function of three parameters: The number of tests perform per unit time, the effectiveness of the test, and the probability of the person that tests positive complying with the quarantine order.
Let's first define all the widgets that we will need.
# Average infective contacts per day for symptomatic
contact_sym = pn.widgets.FloatSlider(
name="infective contact symptomatics (βs)",
start=0,
end=2,
step=0.05,
value=0.1,
)
# Average infective contacts per day for asymptomatic
contact_asym = pn.widgets.FloatSlider(
name="infective contacts asymptomatics (βa)",
start=0,
end=2,
step=0.05,
value=0.2,
)
# Average infective contacts per day for asymptomatic in quarantine
contact_asym_q = pn.widgets.FloatSlider(
name="infective contacts quarantined (βaq)",
start=0,
end=2,
step=0.05,
value=0,
)
# Incubation period
incu_time = pn.widgets.FloatSlider(
name="incubation period (Ɣ⁻¹) [days]", start=0, end=30, step=1, value=5
)
# Fraction of symptomatic patients
p_sym = pn.widgets.FloatSlider(
name="symptomatic fraction", start=0, end=1, step=0.01, value=0.14
)
# Rate of testing (fraction of population per day)
k_test = pn.widgets.FloatSlider(
name="rate of testing (fraction of population per day)",
start=0,
end=1,
step=0.00001,
value=0.0001,
)
# test efectiveness
delta_eff = pn.widgets.FloatSlider(
name="test effectiveness (δeff)", start=0, end=1, step=0.05, value=0.70,
)
# test efectiveness
p_comp = pn.widgets.FloatSlider(
name="probability of quarentine compliance",
start=0,
end=1,
step=0.05,
value=0.80,
)
# Time to recovery
recovery_time = pn.widgets.FloatSlider(
name="time to recovery (µ⁻¹) [days]", start=7, end=30, step=1, value=15,
)
# Population size
log_pop = pn.widgets.FloatSlider(
name="log pop. size", start=1, end=10, step=0.1, value=np.log10(3E8)
)
# Integration time
int_time = pn.widgets.IntSlider(
name="integration time [days]", start=100, end=1000, step=1, value=100
)
# Plot cumulative infected people
infect_toggle = pn.widgets.Toggle(
name='total infected', button_type='success'
)
Having defined these widgets we will now define a function that uses the values of these widgets, computes the dynamics, and returns an interactive plot of each of the populations.
@pn.depends(
contact_sym.param.value,
contact_asym.param.value,
contact_asym_q.param.value,
incu_time.param.value,
p_sym.param.value,
k_test.param.value,
delta_eff.param.value,
p_comp.param.value,
recovery_time.param.value,
log_pop.param.value,
int_time.param.value,
infect_toggle.param.value,
)
def SEIIIR_model(
contact_sym,
contact_asym,
contact_asym_q,
incu_time,
p_sym,
k_test,
delta_eff,
p_comp,
recovery_time,
log_pop,
int_time,
infect_toggle,
):
"""
test
"""
# Define time (in days) for integration
t = np.linspace(0, int_time, 500)
# Set parameters
pop_size = 10**log_pop
beta_s = contact_sym / pop_size
beta_a = contact_asym / pop_size
beta_aq = contact_asym_q / pop_size
gamma = 1 / incu_time
gamma_s = gamma * p_sym
gamma_a = gamma * (1 - p_sym)
theta = k_test * delta_eff * p_comp
mu = 1 / recovery_time
# Pack parameters into tuple
args = (
beta_s,
beta_a,
beta_aq,
gamma_s,
gamma_a,
theta,
mu,
)
# Set initial conditions for population
init_infection = 1e-2
x0 = np.array(
[
1 - init_infection, # S
0, # E
init_infection / 3, # I_s
init_infection / 3, # I_a
init_infection / 3, # I_aq
0, # R
]
) * pop_size
# Numerically integrate system
x = sp.integrate.odeint(system_rhs, x0, t, args=args) / pop_size
# Bokeh plot
# Set colors for curves
colors = bokeh.palettes.d3["Category10"][10]
# Split populations
S, E, I_s, I_a, I_aq, R = x.transpose()
# Setup plot
p = bokeh.plotting.figure(
plot_width=500,
plot_height=300,
x_axis_label="time (days)",
y_axis_label="fraction of population",
)
# Populate plot with glyphs
p.line(t, S, line_width=2, color=colors[0], legend_label="S")
p.line(t, E, line_width=2, color=colors[1], legend_label="E")
if infect_toggle:
p.line(
t, I_s + I_a + I_aq,
line_width=2,
color=colors[6],
legend_label="I"
)
else:
p.line(t, I_s, line_width=2, color=colors[2], legend_label="I_s")
p.line(t, I_a, line_width=2, color=colors[3], legend_label="I_a")
p.line(t, I_aq, line_width=2, color=colors[4], legend_label="I_aq")
p.line(t, R, line_width=2, color=colors[5], legend_label="R")
# Place legend
p.legend.location = "center_right"
# Apply PBoC plot format
covid19.viz.pboc_single(p)
return(p)
Let's now arrange all parameter widgets into a dashboard to start exploring the parameters space.
pn.Column(
pn.WidgetBox(
pn.Row(
pn.Column(contact_sym, contact_asym, contact_asym_q),
pn.Column(incu_time, p_sym, recovery_time),
pn.Column(k_test, delta_eff, p_comp)
)
),
pn.Row(
SEIIIR_model, pn.WidgetBox(pn.Column(log_pop, int_time, infect_toggle))
)
)
An interesting quantity to explore under different scenarios what would be the maximum number of infected people at any point. In particular this is an interesting quantity to explore as a function of the number of tests available for the population. Let's again set an interactive plot, recycling some of the widgets generated for the previous plot.
# Define range of testing rates (fraction of pop. per day)
k_test_range = np.linspace(0, 1, 100)
# Establish function dependencies
@pn.depends(
contact_sym.param.value,
contact_asym.param.value,
contact_asym_q.param.value,
incu_time.param.value,
p_sym.param.value,
delta_eff.param.value,
p_comp.param.value,
recovery_time.param.value,
log_pop.param.value,
)
def max_infection(
contact_sym,
contact_asym,
contact_asym_q,
incu_time,
p_sym,
delta_eff,
p_comp,
recovery_time,
log_pop,
):
# Define time (in days) for integration
t = np.linspace(0, 500, 1000)
# Set parameters
pop_size = 10 ** log_pop
beta_s = contact_sym / pop_size
beta_a = contact_asym / pop_size
beta_aq = contact_asym_q / pop_size
gamma = 1 / incu_time
gamma_s = gamma * p_sym
gamma_a = gamma * (1 - p_sym)
mu = 1 / recovery_time
# Initialize array to save maximum number of infected cases
I_max = np.zeros_like(k_test_range)
# Loop through k_test values
for i, k in enumerate(k_test_range):
# Compute theta
theta = k * delta_eff * p_comp
# Pack parameters into tuple
args = (
beta_s,
beta_a,
beta_aq,
gamma_s,
gamma_a,
theta,
mu,
)
# Set initial conditions for population
init_infection = 1e-2
x0 = (
np.array(
[
1 - init_infection, # S
0, # E
init_infection / 3, # I_s
init_infection / 3, # I_a
init_infection / 3, # I_aq
0, # R
]
)
* pop_size
)
# Numerically integrate system
x = sp.integrate.odeint(system_rhs, x0, t, args=args) / pop_size
# Split populations
S, E, I_s, I_a, I_aq, R = x.transpose()
# Add number of infected and save it
I_max[i] = np.max(I_s + I_a + I_aq)
# Set Bokeh plot
p = bokeh.plotting.figure(
plot_width=500,
plot_height=300,
x_axis_label="fraction of population tested per day",
y_axis_label="max. number of infections",
)
# Populate plot with glyphs
p.line(
k_test_range, I_max, line_width=2,
)
# Apply PBoC plot format
covid19.viz.pboc_single(p)
return(p)
As before, let's arrange the widgets into a dashboard.
pn.Column(
pn.WidgetBox(
pn.Row(
pn.Column(contact_sym, contact_asym, contact_asym_q),
pn.Column(incu_time, p_sym, recovery_time),
pn.Column(delta_eff, p_comp, log_pop)
)
),
max_infection
)
Incubation period